<?xml version="1.0" encoding="iso-8859-1"?>
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
"http://www.w3.org/TR/html4/loose.dtd">
<html>
<head>
<link rel="stylesheet" type="text/css" href="newstyle.css">
<meta http-equiv="Content-Type" content="text/html;
  charset=iso-8859-1">
<style type="text/css">
<!--
kbd {
  font-family: "Nimbus Mono L", "Courier New", Courier, Monospace;
}
-->
</style>
<title>Root Rational Fraction program</title>
</head>
<body>

<div id="top">
<h2>Anthony Stone's root-rational-fraction program package</h2>
</div>
<div id="trail">
<p><a href="http://www-stone.ch.cam.ac.uk/">Home</a>
 -> <a href="http://www-stone.ch.cam.ac.uk/programs.html">Programs</a>
 -> <a href="index.html">RRF</a>
 -> Subroutines and functions</p>
</div>

<!-- Index starts here -->
<div id="index">
<p class="indexentry"><a href="http://www-stone.ch.cam.ac.uk/index.html">Home</a></p>
<p class="indexentry"><a href="http://www-stone.ch.cam.ac.uk/personal/personal.html">Personal information</a></p>
<p class="indexentry"><a href="http://www-stone.ch.cam.ac.uk/programs.html">Computer programs</a></p>
<p class="indexentry"><a
href="http://www.ch.cam.ac.uk/staff/ajs.html">Research details</a></p>
<p class="indexentry"><a href="http://www-stone.ch.cam.ac.uk/personal/timf.html">The Theory of Intermolecular Forces</a></p>
<p class="indexentry"><a href="http://www-stone.ch.cam.ac.uk/personal/bibliography.html">Publications</a></p>
<p class="indexentry"><a href="http://www-stone.ch.cam.ac.uk/personal/supplementary.html">
Supplementary material for published papers</a></p>
</div>
<!-- Index ends here -->

<div id="content">

<p>This set of Fortran 90 routines performs arithmetic on numbers
whose squares are rational fractions (possibly negative, so that the
numbers may be pure imaginary). In particular, the theory of angular
momentum yields various numbers of this form (the Wigner 3j, 6j and 9j
coefficients) and routines are provided for calculating them. The
source code is available as a <a
href="http://www-stone.ch.cam.ac.uk/pub/rrf-4.0.tgz">gzipped
tarfile</a> or  <a
href="http://www-stone.ch.cam.ac.uk/pub/rrf-4.0.zip">zip file</a>
(about 25 kbytes).</p>

<p>There is a program called <a href="rrfcalc.html">rrfcalc</a>, included
in the distribution, that uses these 
routines and that can be used interactively to calculate Wigner 3j, 6j
and 9j coefficients and to perform simple arithmetic on them in the manner
of a <a href="rrfcalc.html">desk calculator</a>.</p>

<p>Each RRF (root-rational-fraction) is referred to via a Fortran
variable of <kbd>type(rrf)</kbd>. Such a variable contains a pointer to a
<kbd>type(node)</kbd> variable, which is the head of a circular list in which the
RRF is held in power-of-prime form. It also contains an integer factor
which may hold an intermediate number that has not yet been factorized
into a product of primes. Such intermediate numbers can get extremely
large in complicated manipulations, so it is recommended to use the
gfortran compiler, which provides 16-byte integers (up to 38 decimal
digits). Other compilers, limited to 8-byte integers (18 decimal
digits) will however be adequate for normal use. The program checks
for integer overflow and should stop rather than give erroneous results.</p>

<p>All arithmetic and i/o manipulations are carried out by special
routines provided in the rrf module, and the list elements are not
available directly to a program that uses <kbd>type(rrf)</kbd> 
variables.</p>

<p>Further details of the program, in its original Fortran 77 form,
may be found in the paper by A. J. Stone and C. P. Wood, <em>Computer
Physics Communications</em> (1980) <strong>21</strong>, 213. 
The current version of the program comprises a Fortran module
rrf_module.F90 containing all the basic routines, a module wigner.F90
  containing routines to calculate 3j, 6j and 9j coefficients, the
  calculator program rrfcalc.F90, and programs test3j.f90, test6j.f90
  and test9j.f90 which carry out various tests of the wigner routines.
</p>

<h2><a id="ix">Index</a> of subroutines</h2>
<h3><a href="#initialization">Initialization</a></h3>
<ul>
<li><a href="#start">INIT_RRF</a> </li>
</ul>
<h3><a href="#housekeeping">Housekeeping</a></h3>
<ul>
<li><a href="#clear">CLEAR</a> </li>
<li><a href="#copy">COPY</a> </li>
<li><a href="#rename">RENAME</a> </li>
<li><a href="#exch">EXCH</a> </li>
<li><a href="#unit">UNIT</a> </li>
<li><a href="#iszero">ISZERO</a> </li>
<li><a href="#nonzero">NONZERO</a> </li>
<li><a href="#chsign">CHSIGN</a> </li>
<li><a href="#conjg">CONJUGATE</a> </li>
</ul>
<h3><a href="#arithmetic">Arithmetic and Conversion</a></h3>
<ul>
<li></li>
<li><a href="#mult">MULT</a> </li>
<li><a href="#multi">MULTI</a> </li>
<li><a href="#setf">SETF</a> </li>
<li><a href="#multf">MULTF</a> </li>
<li><a href="#power">POWER</a> </li>
<li><a href="#add">ADD</a> </li>
<li><a href="#accum">ACCUM</a> </li>
<li><a href="#subtr">SUBTR</a> </li>
<li><a href="#equal">EQUAL</a>  </li>
<li><a href="#root">ROOT</a>  </li>
<li><a href="#fromi">INT_TO_RRF</a> </li>
<li><a href="#fromch">PP_TO_RRF</a> </li>
<li><a href="#display">DISPLAY</a> </li>
<li><a href="#to4i">RRF_TO_4I</a> </li>
<li><a href="#char4i">CHAR4I</a> </li>
<li><a href="#write4i">WRITE4I</a> </li>
<li><a href="#toreal">RRF_TO_REAL</a> </li>
<li><a href="#tochar">RRF_TO_CHAR</a> </li>
<li><a href="#list">LIST</a>  </li>
</ul>
<h3><a href="#angular">Angular Momentum Functions</a></h3>
<ul>
<li><a href="#threej">THREEJ</a> </li>
<li><a href="#sixj">SIXJ</a> </li>
<li><a href="#ninej">NINEJ</a> </li>
<li><a href="#three0">THREE0</a> </li>
<li><a href="#regge3">REGGE3</a>  </li>
</ul>
<h3><a href="#errors">Error handling</a></h3>
<h3><a href="#notes">Notes</a></h3>
<hr>

<h2><a id="Initialization">Initialization</a> routine:</h2>
<dl>
<dt>SUBROUTINE <a id="start">INIT_RRF([quiet[,primes]])</a></dt>
<dd>This must be called before any of the other routines. It initialises the
list and sets up a table of primes. QUIET is an optional logical
argument, which if true suppresses the initial program banner.
PRIMES is an optional integer argument, which if present specifies the
number of primes to be generated in the table of primes. Default 200000.
</dd>
</dl>

<h2><a id="housekeeping">Housekeeping</a> routines:</h2>

<dl>
<dt>SUBROUTINE <a id="clear">CLEAR</a>(K) </dt>
<dd>Deallocate RRF K. This should always be used to
release the space occupied by temporary variables before exit from
subroutines -- otherwise the list space will rapidly become full.
There is no automatic garbage collection. </dd>

<dt>SUBROUTINE <a id="copy">COPY</a>(K1,K2) </dt>
<dd>Set RRF K2 equal to K1. (Note that this is not the same as the Fortran
statement K2=K1, which must not be used.) </dd>

<dt>SUBROUTINE <a id="rename">RENAME</a>(K1,K2) </dt>
<dd> Set RRF K2 equal to K1 and clear K1. </dd>

<dt>SUBROUTINE <a id="exch">EXCH</a>(K1,K2) </dt>
<dd>Exchange the RRFs K1 and K2. Use EXCH or RENAME rather than COPY
where possible. </dd>

<dt>SUBROUTINE <a id="unit">UNIT</a>(K) </dt>
<dd>Set K to the RRF value of unity (not the same as K=1). </dd>

<dt> FUNCTION <a id="iszero">ISZERO</a>(K) </dt>
<dd>Returns the value .TRUE. if K represents zero and .FALSE.
otherwise.</dd> 

<dt> FUNCTION <a id="nonzero">NONZERO</a>(K) </dt>
<dd>Returns the value .TRUE. if K represents a nonzero value and
.FALSE. otherwise.</dd> 

<dt>SUBROUTINE <a id="chsign">CHSIGN</a>(K) </dt>
<dd>Change the sign of K (multiply by -1). </dd>

<dt>SUBROUTINE <a id="conjg">CONJUGATE</a>(K) </dt>
<dd>Replace K by its complex conjugate. </dd>
</dl>
<p>Up to <a href="#ix">index</a> </p>

<h2><a id="arithmetic">Arithmetic</a> routines:</h2>

<dl>
<dt>SUBROUTINE <a id="mult">MULT</a>(K1,K2,M) </dt>
<dd>Multiply K1 by K2**M. K1 and K2 must be different variables. 
Division is achieved by having M&lt;0. </dd>

<dt>SUBROUTINE <a id="multi">MULTI</a>(K1,I,M) </dt>
<dd>Multiply RRF K1 by I**M, I, M integers. </dd>

<dt>SUBROUTINE <a id="setf">SETF</a>(N,K) </dt>
<dd>Set RRF K to the value of factorial(N). The factorial table, which
initially contains only factorial(0) and factorial(1), is extended if
necessary. </dd> 

<dt>SUBROUTINE <a id="multf">MULTF</a>(K,N,M) </dt>
<dd>Multiply RRF K by (factorial N)**M. The factorial table is extended if necessary. </dd>

<dt>SUBROUTINE <a id="power">POWER</a>(K1,M) </dt>
<dd>Replace K1 by K1**M. </dd>

<dt>SUBROUTINE <a id="add">ADD</a>(K1,K2) </dt>
<dd>Add K2 to K1, leaving K2 unchanged. K1 and K2 must be different
variables. An additional restriction, which is less of a limitation
than it might seem, is that the result of the addition must itself be
expressible as an RRF. Thus, for example, an attempt to add sqrt(3) to
sqrt(5) will provoke an error message. See Note 2. </dd>

<dt>SUBROUTINE <a id="subtr">SUBTR</a>(K1,K2) </dt>
<dd>Subtract K2 from K1. </dd>

<dt>FUNCTION <a id="equal">EQUAL</a>(K1,K2) </dt>
<dd>The logical function EQUAL has the value .TRUE. if K1 and K2
represent the same number, .FALSE. otherwise. Note that this is not
the same as K1 .EQ. K2. </dd>

<dt>SUBROUTINE <a id="root">ROOT</a>(K) </dt>
<dd>Replace K by its square root. </dd>
</dl>

<p>Up to <a href="#ix">index</a> </p>

<h2><a id="io">Conversion</a> routines:</h2>

<p>These routines convert RRFs to ordinary Fortran variables, or to values
in an A1 character buffer, and vice versa, or print out a character
representation of an RRF. </p>
<dl>
<dt>SUBROUTINE <a id="fromi">INT_TO_RRF</a>(I,K) </dt>
<dd>Express integer I as RRF K. </dd>

<dt>SUBROUTINE <a id="fromch">PP_TO_RRF</a>(M,N, NP, K) </dt>
<dd>Read an RRF from the first N characters of the A1 character buffer M. The
detailed syntax of the RRF representation is given <a href="ppx.html">here</a>. </dd>

<dt>SUBROUTINE <a id="to4i">RRF_TO_4I</a>(K, I1,I2,I3,I4) </dt>
<dd>Express RRF K as (I1/I2)sqrt(I3/I4). I1 and I3 may be negative. I2 and I4
are always positive, unless integer overflow occurs, in which case zero
is returned in I2, and I1, I3 and I4 contain rubbish. </dd>

<dt>SUBROUTINE <a id="char4i">CHAR4I</a>(K, M,M1,M2) </dt>
<dd>Convert the rrf K to 4-integer form as characters in the buffer M of length
M2, starting at position M1. M1 is updated to point at the next available
position in the buffer. If integer overflow occurs, 12 asterisks are returned. </dd>

<dt>SUBROUTINE <a id="toreal">RRF_TO_REAL</a>(K,A,I) </dt>
<dd>Express RRF K as A*i**I, where A is double precision and the integer I has
the value 0 or 1. </dd>

<dt>SUBROUTINE <a id="tochar">RRF_TO_CHAR</a>(K, M, M1,MAX, NP) </dt>
<dd>Express RRF K in A1 character form in the buffer M(MAX), starting at position
M1. The sign term is output as +, -, +i or -i; the exponents of the first
NP primes follow, and further terms appear as '. prime^exp'. NP may be zero.
A value of zero, +1 or -1 is output as 0, +1 or -1 in the sign position.
An error message is printed if there is not enough space in the buffer.
M1 is updated to point at the next available position in the buffer. </dd>

<!--
<dt>SUBROUTINE <a id="list">LIST</a>(K) </dt>
<dd>This is a debugging aid. It prints the actual list entries for RRF
K on unit 6. If it is called with a negative argument -N, it prints
the last N list entries. Each entry occupies three integer locations;
the first contains a pointer to the next entry, the second contains a
prime, and the third is twice the exponent of that prime. The primes
occur in increasing order, and the last entry points back to the
first, for which the 'prime' is always -1. 'Primes' larger than the
square of the largest entry in the prime table may in fact be
composite, but this will not affect the working of the program. </dd>
-->
</dl>

<p>Up to <a href="#ix">index</a> </p>

<h2><a id="angular">Angular momentum functions (module wigner):</a></h2>
<dl>
<dt>SUBROUTINE <a id="threej">THREEJ</a>(J1,J2,J3, M1,M2,M3, X, K) </dt>
<dd></dd>
<dt>SUBROUTINE <a id="sixj">SIXJ</a>(J1,J2,J3, L1,L2,L3, X, K) </dt>
<dd></dd>
<dt>SUBROUTINE <a id="ninej">NINEJ</a>(A,B,C, D,E,F, G,H,I, X, K) </dt>
<dd>These routines set the appropriate vector coupling coefficient in the RRF
K. X is a Fortran integer which should have the value 1 or 2; each argument
(J1 etc.) is interpreted as J1 or J1/2 according as x="1" or 2. </dd>

<dt>SUBROUTINE <a id="three0">THREE0</a>(J1,J2,J3, K) </dt>
<dd>A more efficient equivalent of SUBROUTINE THREEJ(J1,J2,J3, 0,0,0, 1, K) </dd>

<dt> SUBROUTINE <a id="regge3">REGGE3</a>(JP1,JP2,JP3, JM1,JM2,JM3, K) </dt>
<dd>An alternative way to get a 3-j coefficient. jp1=J1+M1, jm1=J1-M1, etc., so all
arguments are integers even if the quantum numbers are not. Since the 3-j
routine needs J1+M1 etc. anyway, this is actually a more efficient way to
define the coefficient required. </dd>
</dl>
<p>Up to <a href="#ix">index</a> </p>

<h2><a id="errors">Error handling</a></h2>

<p>There is a variable <kbd>hard_errors</kbd> in module rrfdata which
is initially set true. With this setting most errors cause the program
to stop after printing an error message. If it is set false, errors
are soft, and cause the variable <kbd>error</kbd> in module
rrflist to be set true after printing the error message. If another
error occurs while <kbd>error</kbd> is true, the program will
treat it as a hard error even if the variable <kbd>hard_errors</kbd>
is false. If <kbd>error</kbd> is set false after an error, the
program will continue, but results should be treated with suspicion.
</p>

<p>Up to <a href="#ix">index</a> </p>
<h2><a id="notes">Notes</a>:</h2>
<ul>
<li>The program should not be used unless exact results are essential, as is
it inevitably much slower than ordinary arithmetic. In particular, addition
may be very slow indeed &#8212; much slower than multiplication. </li>
<li>The embargo on adding numbers if the result cannot be expressed as an RRF
will seem a very severe restriction. In fact, the program has been used
for several fairly elaborate calculations, involving heavy use of 3j, 6j
and 9j coefficients, without ever coming up against the restriction. I would
be interested to learn of a serious problem in which it is an important
restriction. It could be removed in principle, but would involve a complete
rewrite of the program. </li>
<li>Machine-dependent values are set in rrf_module.F90, and
should be checked before compiling the program. </li>
<li>The angular momentum routines all use factorials. The largest available
value is factorial(200), which should suffice for most purposes. An error
message is printed if this limit is exceeded, in which case it is necessary
to recompile with a larger value for the parameter MAXFCT in rrf_module.F90.</li>
<li>The program has been fairly thoroughly tested and has been used
successfully on a number of problems, but some bugs doubtless remain.
I should appreciate being told about them.  </li>
</ul>
</div>
<!--		Ends Here			-->

<div id="foot">
<!--#config timefmt="%d %B, %Y" -->
<p>This page is maintained by
<a href="mailto:ajs1%40cam.ac.uk">Anthony Stone</a>,
and was last updated on <!--#flastmod file="routines.html" -->. </p>
</div>
<div id="w3valid">
<p>
    <a href="http://validator.w3.org/check?uri=referer"><img
      src="http://www.w3.org/Icons/valid-html401" alt="Valid HTML 4.01 Transitional" height="31" width="88"></a>
  </p>
</div>
</body>
</html>
